Evidence for Quadratic Tidal Tensor Bias from the Halo Bispectrum 
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The relation between the clustering properties of luminous matter in the form of galaxies and the 
underlying dark matter distribution is of fundamental importance for the interpretation of ongoing 
and upcoming galaxy surveys. The so called local bias model, where galaxy density is a function 
of local matter density, is frequently discussed as a means to infer the matter power spectrum 
or correlation function from the measured galaxy correlation. However, gravitational evolution 
generates a term quadratic in the tidal tensor and thus non-local in the Eulerian density field, even 
if this term is absent in the initial conditions (Lagrangian space). Because the term is quadratic, 
it contributes as a loop correction to the power spectrum, so the standard linear bias picture still 
applies on very large scales, however, it contributes at leading order to the bispectrum for which it is 
significant on all scales. Such a term could also be present in Lagrangian space if halo formation were 
influenced by the tidal field. We measure the corresponding coupling strengths from the matter- 
matter-halo bispectrum in numerical simulations and find a non-vanishing coefficient for the tidal 
tensor term. We find no scale dependence of the inferred bias parameters up to k ~ 0.1 /iMpc -1 
j ^ and that the tidal effect is increasing with halo mass. While the local Lagrangian bias picture is a 

better description of our results than the local Eulerian bias picture, our results suggest that there 
might be a tidal tensor bias already in the initial conditions. We also find that the coefficients of the 
(N quadratic density term deviate quite strongly from the theoretical predictions based on the spherical 

collapse model and a universal mass function. Both quadratic density and tidal tensor bias terms 

Omust be included in the modeling of galaxy clustering of current and future surveys if one wants to 
achieve the high precision cosmology promise of these datasets. 
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Large Scale Structure (LSS), the large-scale distribution of matter in the Universe, contains a wealth of information 
about the history and composition of the Universe as well as fundamental physics. For instance, LSS has the potential 
to constrain neutrino masses or modifications of gravity, which however requires percent level accuracy for the theory 
and observations. Besides gravitational lensing, which is sensitive to the total matter distribution, the positions of 
galaxies are the main observable and tool to infer the underlying matter distribution. They have the advantage of 
higher statistical power relative to weak lensing surveys. Ongoing and upcoming LSS surveys such as BOSS, BigBOSS, 
EUCLID, DES will provide an unprecedented quality of galaxy clustering data, which needs to be properly analyzed. 

A crucial step in the analysis of galaxy surveys is to connect the distribution of the tracer to the underlying 
distribution of matter. The first step in this logical chain is the realization that galaxies form preferentially in the 
potential wells of collapsed dark matter halos, where the hot gas can cool sufficiently fast [TJ [2]. This leads to the 
question of how the clustering properties of dark matter halos relate to the clustering of matter in general. The answer 
to this question is usually phrased in terms of a relation between the overdensities in these two fields and is dubbed 
a halo biasing scheme. With the ever increasing computing power it is in principle possible to generate templates for 
the survey analysis for standard ACDM and even modified gravity models using iV-body simulations. This approach 
becomes very expensive when it is to be used in Markov-Chain-Monte-Carlo parameter inference methods and does not 
provide insight into the underlying clustering properties. We thus consider it important to understand the properties 
of halo clustering by testing theoretical prescriptions on simulations with the final goal of devising analytical and thus 
easily evaluable models for the survey analysis. 

The so called local biasing model [3l 0] , where galaxy and halo density is a function of local matter density, has 
been the most popular model used in previous work. In the simplest version one adds another contribution that 
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scales quadratically with density, the quadratic density bias term. Recent work has argued, based on symmetry and 
analyticity arguments, that there are additional terms not included in the local bias model that appear in the power 
spectrum and are formally at the same order as that of quadratic bias [5] . One of these terms is quadratic in density 
but non-local and can be written as the square of the tidal tensor. The first goal of this paper is to provide additional 
theoretical motivation for inclusion of this term in the analysis of galaxy clustering. 

The working horses in LSS analysis are the two-point functions, the correlation function and the power spectrum. 
Purely Gaussian, linear fields are completely characterized by their two point function. However, non-linear phe- 
nomena in galaxy and halo formation as well as non-linear gravitational clustering can generate the full hierarchy of 
n-point functions. These higher order statistics might be difficult to measure in the sky due to non-trivial survey 
windows and redshift space distortions, but they can be easily extracted from ./V-body simulations. The simplest 
statistic beyond the power spectrum is the bispectrum 

(6{k 1 )6(k 2 )5(k 3 )) = B(k ll k 2 ,k 3 ) (2^) 3 <5 (D) (fei + fc 2 + fc 3 ), (1) 

which is well suited for the study of next-to-leading order effects in cosmic density fields [rjHTU] . While these contribute 
only loop terms to the power spectrum they are the leading order terms for the bispectrum, which vanishes for purely 
linear Gaussian fields. 

The second aim of this study is to probe the halo bispectrum for Gaussian initial conditions in the low-fc regime, in 
order to extract the two terms that lead to a quadratic coupling between the number density of collapsed objects and 
the long wavelength matter fluctuations, quadratic density bias and quadratic tidal tensor bias. We present a study 
of their scale and mass dependence using TV-body simulations and we compare the numerical results to the theoretical 
expectations based on simple halo bias models. 

This paper is organized as follows. In Sjlljwe review the standard formulation of the bias model and discuss possible 
extensions. §III| describes the simulations, the bispectrum measurement and data reduction as well as the parameter 
estimation. The results are presented in §IV| Finally, in fjVJwe discuss our findings and their implications as well as 
possible directions for future investigation. 

II. THE BIAS MODEL: LOCAL AND NON-LOCAL FORMS 

A. Standard Formulation: Local Bias 

The formation of galaxies and their host dark matter halos is a complicated highly non-perturbative process. It 
is, however, reasonable to assume that certain properties of collapsed objects, for instance their number density, are 
related to the coarse grained underlying matter density field in the same region of space. Neglecting complications 
arising from gas physics, the number density of collapsed objects can be written as a functional of the underlying 
matter density perturbation [3] 

6 h (x,r))=F[6(x', V )]. (2) 

On large scales, this functional is commonly approximated by a local and linear bias model <5h (x, rj) — b\5 (x, rf) (plus 
generally a noise term which we will avoid in this paper by only looking at cross-correlations with mass) . The next step 
is to give up on linearity and to introduce the second order local bias model <5h {x, rf) = b\8 (x, rj) + &2<5 2 (x, rj). This 
model has been well studied in the literature in combination with Standard Perturbation Theory (SPT) and leads to 
non-trivial renormalizations of the leading order bias parameter . Measurements of the quadratic bias parameters 
of this model in the bispectrum have lead to contradictory results, which raises doubts about the completeness of the 
model [T2] (see also [13]). 

B. Tidal Terms 

The power series expansion of the functional presented above is certainly overly simplified and one should consider 
whether other terms could influence the number density of collapsed objects. As proposed in [5], the environmental 
dependence of halo formation could lead to a dependence on the tidal field, as quantified by the tidal tensor 

Sij{x,rj) = didj$(x,r}) - ^ K) <5(a:, 77). (3) 

Note that we absorbed the constants in the Poisson equation into the gravitational potential 'V 2 Q(x,r)) = 6(x,rj) 
and subtract out the trace from the tidal tensor because it is degenerate with the density field. The corresponding 
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expression in Fourier space is given by 

Sii(k,V)=(^-l4?)6(k,ri)- (4) 

The halo overdensity is a scalar quantity and can thus only depend on scalars. The simplest scalar that can be 
constructed from the tidal tensor is given by s 2 (a;) = Sjj {x)sij (x) which in Fourier space is expressed by the convolution 



(k, n) = J S 2 {k', k - k')S(k', r,)5{k - fe', 77) , (5) 



where we have implicitly defined the kernel 

Sa(gi,ga) = — o o- ( 6 ) 

yiy2 

Following [5J, the halo density field up to second order can be written as 

S h {x, n) = btS(x, v) + b 2 [S 2 (x, rj) - (5 2 (x, 77))] + b s 2 [s 2 (x, rj) - (s 2 (x, , (7) 

where we absorbed prefactors of 1/2 into the bias parameters. We truncated the series at second order, since higher 
order terms influence the bispectrum only through loop corrections, which are believed to be subdominant on large 
scales. As can be easily verified from the above definitions one has (s 2 (ie,77)) = 2/3 (y^5 2 {x, 77)). 



C. Lagrangian Bias 

In the usual Lagrangian bias picture, the galaxy formation sites are identified in the primordial density field, and it 
is assumed that the primordial halo density field at initial time rji can be written as a power series in the primordial 
matter fluctuations. For calculational convenience, the expansion can be rewritten in terms of the linearly extrapolated 
density field S(q,r]) = D(rj)/D(r] i )S(q) 

- E ^r 1 [<%) - <<%)>] - E ^r 1 [*W - <*W)>] , (8) 
1 1 

where b[ L \rj) = (D(r)i) / D(rf)) bf^fa) and q is the Lagrangian position. Here, we introduced the conformal time 
adr/ — dt and the linear growth factor D(r]), normalized to unity at present time. We start the sum from / = 1 
because the I = term vanishes by requiring that the halo field has a vanishing mean. In contrast to this, the Eulerian 
bias model expands the halo density field at a certain point in time in the non-linear matter density field at the same 
time. Motivation for a Lagrangian nature of halo bias comes from the peak model [T4lI17] . where the peaks of the 
primordial density field are associated with the formation sites of protohaloes. 

It now remains to connect the Lagrangian density fields to the observable Eulerian ones. The continuity equation 
for halos requires 

[1 + 5 h (x, 77)] d 3 x = [1 + 5 h (q)} d 3 q, (9) 

where 

[l + 6(x,r))} d 3 x = d 3 q (10) 

is the continuity equation for the underlying dark matter field. Note that the Lagrangian density field we work with 
is always the linear density field, because this was the definition of the bias expansion in Eq. Q - more general 
expansions are possible in principle. 

The Lagrangian and Eulerian positions are related by x(q, rj) = q + ^f(q, rj), thus up to third order in the density 
field we have 

*(g) = S{x,n) - *i{q,ri)di6(q) + \^i{q,v)^i(9,v)W{Q)- ( U ) 
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In contrast to Eulerian perturbation theory, where the density is the central quantity, Lagrangian Perturbation 
Theory (LPT) has the displacement field \& as the central dynamic quantity and the density fields are only derived 
quantities (for the basics of LPT and its relation to SPT see Appendix [A]). For simplicity, we will focus on the matter- 
only Einstcin-dc Sitter Universe for the theoretical calculations throughout this paper. Using the Lagrangian bias 
expansion j8]) in (T9J) and expressing the Lagrangian position field in terms of the Eulerian position, we have (see also 



4^,77) = (1 + 6^(77)) W8(x,r,) 



-b^( V )S 2 (x,rj) + \b%\r,)^8 2 { X ,rj) - (^5 2 (x, 77))] - b^\vM9,v) • VS(x, V ) + {2) 5(x, v ). 



(12) 



Reorganizing the terms in order to have the first order bias multiply the full second order matter density field, we 
obtain 



S h (x,ri) = (l + &!%)) r?) + {2) S(x, v )) 



(L) 



(13) 



where we used that the second order mass density (in SPT and LPT) can be written in configuration space as (see 
e.g. 
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(2) <5(a:,77) = — {1) S 2 (x,i]) - ^{x,rj) ■ V5(x,rj) + -s 2 (x,i]). 



(14) 



We see that the functional form of the above result agrees with Eq. but complements it with a dynamical 
perspective. In particular, we see that even in the absence of a tidal tensor bias in the initial conditions, such a term 
will be generated by subsequent gravitational evolution. While our result in Eq. ( 13 1 is consistent with the relations 
presented in |18j , we have reorganized terms in order to have the first order bias multiplying the second order density 
field Eq. (14) and to make the formal equivalence with the phenomenological picture presented above more obvious. 
Furthermore, there are no b^ factors, because we fixed the original Taylor series for 8h(q) to have vanishing mean. 



D. Coevolution of Halos and Dark Matter 

In this subsection we will consider an (at first glance) Eulerian approach to halo clustering and consider the 
coevolution of the coupled halo-dark matter fluid (For a similar approach in combination with resummed perturbation 
theory see |20|). Assuming vanishing velocity bias = v and conservation of halo number, we can write down a 
coupled system of differential equations for the matter and halo fluid [38], namely the continuity equation for halos, 
the continuity equation for the matter and the combined Euler and Poisson equations for matter 

s' h (k, v ) + e(k, v) = -J 44 a(k', k - k')e(k', v )s h (k - k', v ) (15) 



(2 



S\k, 77) + e(k, V) = ~J a(k', k - k')6(k', V )S(k - k', rf) 



(16) 



9'(k, 77) + U(r,)6(k, 77) + ln m (r,)H 2 (v)S(k, »/)=-/ * - k')0{k', V )6(k - k', rf) (17) 



2 mV " v,y-v-.-,/ j (2?r 

Here we introduced the velocity divergence 9(x) = V • v(x). The matter equations at first order are solved by 
^6{k,rf) = — H(r])^S(k, 77) and ^S(k,rf) = D(rf)^'So(k), where ^So(k) is the present day linear overdenisty. The 
second order solution for the matter yields 

(2 



(2) 



S(k )V )= J |^^(fc',fc-fc') (1) 5(fc',77)W<5(fc-fc',77) (18) 
6{k,rf)=-U{rf) J ^G 2 (fc',fc-fc') (1) '5(fe','7) (1) '5(fc-fe',77) (19) 



5 



with 



where the second order SPT mode coupling kernels are defined as [3T] 

5 2 

*2 (91,9a) =^a(qi,q 2 ) + -/3(gi,q 2 ) (20) 
3 4 

G 2 (91,92) =^a(9i,92j + =/3(9i,92j (21) 
/ s (91 + 92) • gi al n 1 / , ^2 91 • 92 /orA 

"(9i,92 = 2 /H9i,92 = o (9i +92) 2 2 (22) 

91 2 ^ 

The first order equation for the halos can be solved using the local bias ansatz ^8h.(k, rj) = b^\rj)^5(k, 77), which 
then gives the time evolution of the first order bias as 

&!%)-! _ DM (o* 

This relation known as debiasing, i.e., at late times the bias converges to unity and halo and matter density field 
agree [22] . 



Solving Eq. ( 15 ) at second order using the second order matter solutions, we obtain 

^S h (k lV ) = ^6 h (k, Vi ) + b^( V ) J ^F 2 (k',k-ki)5(k', v )5(k-k>, v ) 



4 ft™. , .\ f d 3 k' 



(2tt) s 

21 (b<?\ V )-l) j 1^6(k>,r,)8(k-k',r,) (24) 
f (ftf'W-l) / ^S 2 (k>,k-k')5(k',ri)5(k-k',Ti), 



where we assumed D(j]i) <C D(rj). Here, we have isolated the part proportional to the second order matter field in 
the first line. Thus, the dynamical evolution naturally introduces a 6 2 (x) and a s 2 (x) term, even in the absence 
thereof at some initial time rj{. Translating back to real space, we see that this has the same functional form as 
Eq. Q and agrees with the Lagrangian bias picture. The equivalence is even more obvious if we have ^S^x, rj\) = 

62 (r)i)^S 2 (x, r]i)/2 = b < ^\r/)^S 2 (x, rj)/2 in correspondence to the Lagrangian bias picture discussed above (see next 
subsection for a relation between the parameters of the models discussed here). We note in retrospect that, while the 
calculation took an Eulerian form, the specification that galaxies formed at an early time tracing the initial density 
field (or at least the location where they will form is determined this way), and then just follow gravity is actually 
the same one made in the Lagrangian calculation. 



E. Relation between Eulerian and Lagrangian bias parameters 

As we have seen above, the Lagrangian bias model, the coevolution model and the educated guess of [5] all lead to 
the same functional form for the halo density field if one identifies the parameters of the model as 

&! = &<%) =1 + 6^(77) (25) 

b2 = \bf\ n )=^bf\ri)+ 1 -bf\r 1 ) . (26) 

This identification agrees with the one of the spherical collapse picture [53]. Note however that our results are not 
relying on the spherical collapse dynamics. For the prefactor of the tidal field scalar we have 

^=-^(^^-1)=-^^. (27) 

In this model there is no tidal field bias in Lagrangian space, consistent with the spherical collapse model, although 
in the ellipsoidal collapse model [24Tf27] such a term would be allowed. 

When comparing our measurements to theoretical bias models, we consider the bias derived from the Sheth-Tormen 
[28j mass function with parameters p = 0.15 and q = 0.75 that were shown to be in good agreement with first order 
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bias in iV-body simulations [25], although we checked that the predictions do not differ much from those using the 
original values in [2 8) . The Lagrangian bias parameters are then given by the first and second derivative of the mass 
function n(M) with respect to a long wavelength background fluctuation Si. For universal mass functions these can 
be rewritten as derivatives with respect to the peak height v(M,r]) = 5 2 /a 2 (M,rj) 



n ao\ n o c av 



2 n OS? n 8}. dv 2 n S 2 dv ' 



1 d 2 n 4 v 2 d 2 n 2 v dn 

where 5 C — 1.686 is the critical density for spherical collapse. The derivatives of the Sheth-Tormen mass function read 

1 dn qv — 1 p 



ndv 1v v{l + {qv)P) 

1 d 2 n p 2 + vpq {qv) 2 — 2qv — 1 
ndv 2 ~ v 2 (1 + {qv)P) + 4^2 ' 



(30) 
(31) 



The mass dependence of the bias function is presented in Figure H] and discussed in £ IV 



F. Bispectra 



The leading order contribution to the bispectrum arises from quadratic terms in the fields. Higher order couplings 
enter only through loop corrections, which gain importance for high k. This is equivalent to the situation in the power 
spectrum, where linear terms in the field are dominant on large scales and loop corrections from quadratic and higher 
order terms gain importance for high k. 

The tree-level matter bispectrum in SPT is given by 



Bmn,m(fci, fe, As 3 ) = 2P(k 1 )P(k 2 )F 2 (k u k 2 ) + 2 eye. 



(32) 



where eye. symbolizes the two cyclic permutations of the fc-vectors in the power spectrum and mode coupling kernel. 
From an observational point of view the halo auto bispectrum -Bhhh is probably the most appealing statistic. Unfor- 
tunately it is suffering from shotnoise, which might deviate from its fiducial Poisson form l/n [3D]. Besides the halo 
auto power spectrum, there are two halo-matter cross bispectra i? m hh and -B mm h, where either one or two matter 
fields are correlated with two or one halo fields, respectively. One further needs to state whether the cross bispectra 
are symmetrized over the k modes or not. 

Our focus is not on observability but on understanding the clustering properties of dark matter halos in A^-body 
simulations and devising a theoretical framework that can later be used to analyze real data. Thus, we will consider 
the un-symmetrized matter-matter-halo cross bispectrum defined as 



B 



(unsym) 
mmh 



(fci,fe 2 ,fe 3 ) (2^) 3 «5( D )(fe 1 + fc 2 + fc 3 ) = (S{ki)S(k 2 )S h (k 3 )) , 



(33) 



where the halo density field is always on the fc3-mode. This particular configuration might not be the one with the 
highest signal to noise ratio but it has a couple of quite useful properties for our study: a) The cross-bispectrum 
does not suffer from a spurious shotnoise contamination and is thus a clean probe of the clustering of halos b) the 
functional simplicity of the second order bias contributions in terms of ki , k 2 and the enclosed angle p where the b 2 
and b s 2 contributions are basically orthogonal (see below) . The latter should allow for a clear distinction between the 
standard second order bias picture or its possible extensions discussed above. 

As noted above, all of the models under consideration share the same functional form Q for the second order halo 
density field. Only the standard quadratic bias model has b s 2 = 0. Therefore, we can write down the bispectrum as 



B 



(unsym) 
mmh 



(fc 1 ,fc 2 ,fe 3 ) - b x B 

mmm 

(k 1 ,k 2 ,k 3 ) = 2P(k 1 )P(k 2 ) 



(34) 



Note that we use a parametrization where the factors of 1/2 are absorbed into the bias parameters for simplicity. We 
will restore these prefactors only at the very end, when we are comparing our bias measurements to the theoretical 
bias functions. A non-vanishing b s 2 in the above equation would be a clear evidence for a dynamical biasing picture. 
After dividing by the two matter power spectra, the remaining quantity is a function of the angle fi only, which 
simplifies the combination of information from several scales k\ and fc 2 . 
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FIG. 1: Matter (black points) and halo-matter-matter (red points) bispectra as a function of triangle shape for configuration 
ki = 0.052 /iMpc -1 &2 = 0.06 /iMpc -1 . The black solid line shows the tree level prediction for the matter bispectrum, the 
red solid line has 6i only and dashed and dashed-dotted lines are adding &2 and b a 2 . The list of bias parameters behind the 
theoretical cross bispectra in the legend indicates the parameters that were considered for the corresponding curve. The error 
bars are estimated from the box-to-box variance of the bispectrum measurement. Note that this is only a small fraction of the 
total bispectrum information that our simulations contain, and the log scale also de-emphasises what are actually significant 
differences between the fit with and without s 2 (these will be highlighted later). 



b% A&i 6 2 A& 2 



Ab s 2 M[fe -1 M S ] 



I 1.142 0.002 

II 1.409 0.002 

III 1.954 0.004 

IV 2.889 0.010 



-0.37 0.01 
-0.42 0.01 
-0.12 0.02 
1.25 0.03 



-0.07 0.03 9.68 x 10 12 

-0.21 0.04 2.90 x 10 13 

-0.38 0.06 8.58 x 10 ls 

-0.63 0.11 2.48 x 10 14 



TABLE I: Best fit bias parameters, their errors and mean mass for our four mass bins. The bias parameters are compared to 
the theoretical bias functions in Figure [4] The first order bias &i is extracted from the halo-matter cross power spectrum and 
the second order bias parameters are inferred from the cross bispectrum. 
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bill IV 
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FIG. 2: Convergence of the measured &2 (upper panels) and b g 2 (lower panels) parameters with increasing maximum fc-mode 
for the four mass bins. The horizontal red and blue lines show the constraints obtained for our fiducial fc max = 0.07 /iMpc -1 . 
The pivot data point is highlighted by the gray shaded region. 



III. SIMULATIONS & BISPECTRUM ESTIMATION 



A. The Simulations 



We are studying the cosmic density field in a suite of 11 dark matter only simulations with box size of L = 
1600 ft. _1 Mpc, which are an extension of the simulations described in [31. The ACDM model is based on best-fit 
parameters inferred from the WMAP 5-year data release [32]. Thus, we adopt a mass density parameter fi m = 0.279, 
a baryon density parameter = 0.0462, a Hubble constant h = 0.7, a spectral index n s — 0.96, and a normalization 
of the curvature perturbations of A 2 ^ = 2.21 x 10~ 9 at the pivot point k = 0.02 Mpc . This normalization leads to a 
present day fluctuation amplitude of a$ ~ 0.81. The initial conditions are set up a redshift z\ — 99. The gravitational 
evolution of the N p = 1024 3 particles is integrated using the publicly available Gadget2 code [33]. Simulations 
size, particle number and matter density parameter yield a particle mass of 3 x 10 11 1%~ 1 Mq. Dark matter halos are 
identified using a Friends-of-Friends halo finder with a linking length of 0.2 times the mean inter particle spacing. Only 
halos exceeding 20 particles are considered for our analysis, corresponding to a minimum halo mass of approximately 
6 x 10 12 h^ 1 M ( ? ) . We consider four mass bins, each spanning a factor of three in mass. For the estimation of the 
statistics, particles are interpolated on a N c — 512 grid using the Cloud-in-Cell (CIC) algorithm and the gridded 
density field is corrected for the window of the grid. The matter-matter-halo and matter bispectrum are measured 
for low /c-modes k < 0.12 /iMpc^ 1 in the simulation output at redshift z = 0. The bispectrum measurement scales as 
the sixth power of the number of grid cells per dimension, which makes it computationally very expensive to extract 
the full bispectrum information at higher k. It would still be interesting to extent the measurement to higher k in 
the future to determine the scale of breakdown of the tree level calculation. 
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bin III 



bin II 




















bin IV 



o.o 

A" 




0.0 



--■ b 2 L (p) 

— bzLoQi) + k,2L. 2 (fi) 

o combined k max — 0.07 /?Mpc~ 

10 05 



FIG. 3: Residual shape dependence of the halo bispectrum for our reduced bispectrum defined in Eq. (351. The blue data 
points with error bars show the result of the combined reduced bispectrum defined in Eq. ( |37[ ) including all the configurations 
up to fc max = 0.07 ftMpc -1 . The horizontal dashed line shows the model including 62 only, the solid blue line shows the model 
including both 62 and b 3 2. 



B. Bispectrum estimation and data reduction 



The bispectrum modes must satisfy the triangle condition k\ + fc 2 + k^ = 0, thus the shape of the bispectrum is 
fully specified by two lengths and one angle, which we choose as &i,fc 2 and fj, = fei • fc 2 /fci£; 2 . Since 8(x) is a real 
valued field, the Fourier modes have to satisfy S*(k) = S(—k). Consequently, the imaginary part of the bispectrum 
cancels when we add S(ki)S(k2)S(k3) and 6(— ki)6(— fc 2 )<5(— k^). We consider bins that are logarithmically spaced in 
k\ and ki and linearly in /i. We add all the bispectrum amplitudes that fall into the bin centered at (fci, fc 2 , /i). 

As a first step in our analysis we subtract out the b\ contribution proportional to the matter bispectrum and divide 
by the power spectrum measured in the same simulation and fc-bins to cancel part of the cosmic variance 

^ ( } = fla m) (fci,ib 2 ,/i)-fcifl mmm (fei,fc 2 ,M) - (35) 

2P mm (A :i )P mm (fc 2 ) 



Comparing to Eq. (34), the resulting statistic should be a function of fi only. The hat is used to mark quantities 
estimated from the simulations. This quantity is distinct from the the usual definition of the reduced bispectrum since 
the power spectra in the denominator do not depend on £3 and can thus be controlled by limiting k\ and fc 2 . The 
first order bias parameter b\ is estimated from the halo-matter cross power spectrum on large scales 0.015 /iMpc -1 < 
k < 0.03 hMpc" 1 , where the linear bias model Ph m = fri-Pmmis believed to be accurate. We show the fc-dependence 
of Phra(k) I P mm (fc) and the inferred bias parameters in Fig. Ul Note that the resulting statistic depends only on the 
magnitude of the ki and fe 2 modes, so that we can ensure the validity of the perturbative expansion by limiting these 
modes accordingly. 
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FIG. 4: Left panel: Mass dependence of the bias parameters and theoretical predictions. The points with error bars are our 
best fit parameters for 61 , 262 and b s 2 . The numerical values of the data poi nts a re given in Table [i] The curves show the 
corresponding theoretical bias functions as calculated using the relations in j jll E[ The measurements for 61 are in a good 
agreement with the theory, there is a clear deviation for the b a 2 and b2 measurement for the two central mass bins. Right panel: 
Ratio of the simulation halo matter and matter power spectra Phm(k) / Pmm{k) and first order bias parameters inferred using 
the data points highlighted by the shaded region. 



When showing the reduced data as a function of /i only, we reduce them as follows 



a _ ^ / M{k l M^)-M{n) \ 

XM = hX ssft^o ) ( } 



A%h (s™) 7 (38) 

for each \x. 

The cosmic variance of the bispectrum estimates could in principle be measured from the standard deviation between 
our simulation boxes. Given the small number of boxes this approach is bound to give a very noisy estimate. Since we 
are using the error for the weighting of the modes, we would like to avoid a spurious upweighting of modes which by 
chance have a low simulation to simulation variance. For this purpose we prefer a smooth error estimate. As shown 
in |34j the variance of the bispectrum is given by 

,2 „ , x V{ 1 „ „. x „ „. , ( n „ , , 1 



AB 1 2 nmh (fc 1 ,fc 2 , j u) = s 123 — -— rP mm (fci)P mm (/c 2 ) P hh (fc 3 ) + — , (39) 
V123 (27r) J V n h J 

where Vf — (2ir) 3 /L 3 is the volume of the fundamental cell, fit is the number density of the tracer and S123 takes on 
values of 1, 2 and 6 for general, isosceles and equilateral triangles, respectively. The norm volume is given by 

V 123 = -^kfk 3 2 ( d In kf d M (40) 

for our bins, which are logarithmically spaced in fc l5 fc 2 and linearly spaced in /i. The above just quantifies the 
diagonals of the covariance matrix between the different triangle shapes and scales, but the correlations between 
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different triangies are believed to be small [31]. When calculating the error on the reduced bispectrum, we focus on 
the error contribution from the matter-matter-halo bispectrum described above and thus have 

AM( fel)fc2 , M )- 2Pinm(fci)Pinni(fc2) (41) 
This procedure could be clearly improved by modeling or measuring the full covariance matrix. 



C. Bias Estimation 



We can now estimate b 2 and b s 2 minimizing 

2 _ ST ST ( M(k 1 ,k 2 ,ii) - b 2 L (p) - b s 2L 2 ([i) \ . . 

x ~kr\ AM (fcl ,, 2 ,,) ) ^ 

where Lq(/i) — 1 and L 2 (/i) — (/i 2 — 1/3) are the zeroth and second order Legendre polynomials, which form an 
orthogonal set on [—1, 1]. The A:-sums are performed over k mm < k\ < fc max , k m - m < k 2 < Kk\ such that (1 — K)ki < 
^3 < (1 + K )ki and we use k = 3/4 for definiteness. 



Defining the cosine 



, \ t / / i n\ V- V- a(ki,k 2 ,ti)b(k 1 ,k 2 ,n) , . 

l (k 1 , k 2> m), b(ki,k 2 ,n)) := ^ ^ (43) 



we obtain for the best fit parameters 



- (L 2 ,L 2 )(L ,M) (L ,L 2 ) (L 2 ,M) .... 
b 2 ^ £ (44) 

- (Lq,L 2 }(L ,M) (L ,Lo)(L 2 ,M) .... 
s2 = A ' ' ^ 

where A = (L a ,L ) (L 2 ,L 2 ) — (L 0l L 2 )' 2 . Note that from x 2 minimization one obtains an equivalent expression for 
the cosine of the reduced data 

:= 2^7=2— I ( 46 ) 
V AM (m) 

where for a(fJ,),b(fi) we have (a(fj.),M(ki,k 2 ,fi)) = (a(/i), Af (//)) and (a(/i). 6(/i)) = (a(fx) , b(fx)) ^. 

The bias parameter &! is estimated from the halo-matter cross power spectrum on large scales/low k, where loop 
corrections are believed to be unimportant. Furthermore there is no shotnoise contamination in the cross power 
spectrum. We use b\ measured from the halo-matter cross-power spectrum and use it for the bispectrum calculation, 
because it has a small statistical error. There is a correlation between the b\ used for the matter bispectrum subtraction 
and the inferred b 2 and b s 2 parameters. The matter bispectrum and the residual considered for our reduced bispectrum 
should in principle not be degenerate, since the matter bispectrum has additional Li(/i) terms and since the coefficients 
of the three Legendre polynomials in the matter bispectrum depend on the ratio k 2 /k\. However, applying our fitting 
method to the matter bispectrum we obtain b 2 and b s 2 constraints of order unity. Thus an error of A&i in the first 
order bias leads to a shift of roughly — Abi for both b 2 and b s 2. The statistical errors on the second order bias 
coefficients are roughly a factor of five larger than the statistical error on the corresponding first order bias, such that 
this effect is not very important from the statistical point of view. We would only need to worry about it if the first 
order bias coefficients in the power spectrum and bispectrum differ. While this may at first glance seem impossible 
(and it probably is), the higher order loop terms renormalize b\ |llj and to our knowledge it has not been explicitly 
shown that the renormalization works the same way in bispectrum as in power spectrum. In this paper we will assume 
that the two are the same. 

The errors on the estimated parameters are calculated from the expected deviation in the x 2 . For a n-component 
parameter vector a the \ 2 i n the vicinity of the best fit parameter set a can be described as 

X 2 (a) = X 2 (a) + \ - ^^I^K- - a d ), (47) 
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where A% 2 = |x 2 ( a ) — X 2 (^)| — 1 corresponds to one sigma errors on the parameters. Thus 



Ab 2 = y/l/(L ,L ) Ab s 2 = ^\/{L 2 ,L 2 ) (48) 

We compared these errors with the standard deviation of the bias constraints obtained from the single realizations 
and found good agreement. 

IV. RESULTS 

Figure [I] shows the matter bispectrum and matter-matter-halo bispectra for our four mass bins for one configuration 
of k\ = 0.052 /iMpc" 1 , k 2 — 0.06 ftMpc -1 as a function of the opening angle fj,. The matter bispectrum is in quite 



good agreement with the theoretical prediction in Eq. ( 32 ) . There is some small tension for positive /x, but this should 
not be a problem for our study since our bias constraints are not relying on this theoretical modeling since we are 
subtracting out the matter bispectrum measured from the simulations. 

We also plot the matter-matter- halo cross bispectrum and the theoretical model of Eq. To visualize the effect 
of the contributions of the bias parameters, we plot the model with vanishing b 2 and b s 2 , with vanishing b s 2 and 
the full model. The non-vanishing bias parameters for the theory lines were chosen according to our best fit model 
discussed below. While the fc-values were chosen quite high to reduce the errors, there is still too much scatter in the 
data points to decide whether the model with and without b s 2 gives a better description of the data. This motivated 
the careful combination of all the information available by weighting the modes accordingly, as described above. 

Our main reason to study lowest order non-linear biasing in the bispectrum was that these terms are the leading 
order terms on large scales. As one includes higher momenta, loop corrections gain importance and need to be modeled 
accordingly. To asses the importance of higher order corrections and to show the convergence of our fitting procedure, 
we perform our parameter estimation as a function of the maximum wavenumber and show the results in Figure [2] 
The error bars clearly shrink as we go to higher k and the inferred bias parameters are almost always consistent with 
our fiducial result shown by the horizontal lines. In these plots we also show the scale up to which we have a complete 
measurement of all the modes by the vertical dashed line and highlight the fc max that we use for the primary reported 
parameter values by the vertical shaded region. 

Figure [3] shows our reduced bispectrum M{fi) defined in Eq. (37) as a function of opening angle /j, and the combined 
errors according to Eq. (38). We overplot the b 2 only model and the model including both b 2 and b s 2. These plots 
show clear evidence for the presence of the tidal term except for the lowest mass bin, for which the b 2 only model 
gives an acceptable description of the angular dependence. 

In Table [I] we give the best fit values of the first and second order bias parameters for our four halo mass bins 
obtained considering all the modes up to fc max = 0.07 /iMpc -1 . Figure [4] shows the bias parameters as a function 
of mass together with the corresponding predictions of the spherical collapse model. Note that we are plotting 2b 2 , 



which corresponds the second order Eulerian bias (see Eq. (26)). The 61 measurements from the halo-matter cross 
power spectrum are in good agreement with the theoretical predictions. The measured b s 2 are slightly lower than the 
theoretical predictions for the two central mass bins, but the trend with mass is well reproduced by the theory. The 
b 2 measurements are less well reproduced by the theoretical bias function, especially mass bins II and III are well 
below the theory. The theoretical predictions for bi and b s 2 are given by first derivatives of the mass function and do 
not reproduce the data perfectly. Thus one would naturally expect some corrections for the second derivatives. The 
disagreement could also be an indication for a failure of the spherical collapse picture at second order. Specifically, the 
fact that the predicted b s 2 disagrees with the measurements suggests that Langrangian b^ is not zero, as predicted 
by the ellipsoidal collapse model [2"4Tf2"T] . 

The right panel of Fig. [4] shows the ratio of halo-matter and matter-matter power spectra used for the inference of 
the first order bias. We are fitting for b\ on large scales to avoid the regime where the non- linear corrections become 
important. These corrections are affecting the highest mass bin quite strongly already starting at k s=s 0.05 ftMpc - . 
The corrections are stronger for higher mass objects, in accordance with the general mass dependence of the second 
order bias parameters derived here, but a full discussion of all the terms entering at one loop level must also include 
the third order terms, which is beyond the scope of this paper. 

As in the power spectrum analysis, the bispectrum is increasingly affected by loop corrections as one increases the 
maximum momentum in the problem. The relevant quantity here is the largest external momentum involved and 
thus one should make sure that all the k's are in the perturbative regime and can be well described by the order of 
perturbation theory considered. Our study focuses on the large scale bispectrum, where a tree level treatment should 
be sufficiently accurate and we have shown that the constraints are both stable and consistent as we increase fc max up 
to 0.1 WVIpc -1 . 
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V. DISCUSSION & OUTLOOK 

In this paper we presented a dynamical motivation for the second order tidal tensor bias term proposed in [5], 
showing that it is naturally generated by gravity even if absent in the initial conditions. We performed a measurement 
of this bias for dark matter halos in simulations, showing clear evidence for the tidal tensor bias, increasing with the 
halo mass. Our results are consistent with the picture in which a significant part of the tidal tensor bias is generated 
by the gravitational evolution, but we also find some evidence that it is present already in the initial conditions. While 
the functional form of the additional terms agrees with the discussion in [5J the dynamical derivation can supplement 
it by a prediction of the time dependence and gives at least qualitative understanding about the mass dependence of 
the bias parameters. 

Our analysis also includes second order density bias bi , which we find to disagree with the theoretical predictions at 
a quantitative level, even if it qualitatively follows the predictions that it should be negative at low mass and increase 
with halo mass. The disagreement calls for a reinvestigation and improvement of the theoretical bias function or 
measurement of the second order bias parameters in the initial conditions, where the Lagrangian bias parameters are 
postulated to describe the density field of the protohaloes. The deviation between the theoretical and measured bi 
parameters does not necessarily mean that the Lagrangian picture is wrong, but rather that the bias functions derived 
from the spherical collapse model and the mass function might not be sufficiently accurate. 

Because these two terms are quadratic, they contribute a loop correction to the power spectrum, so the standard 
linear bias picture still applies on very large scales, even if they contribute at the leading order to the bispectrum. 
If we want to understand halo power spectrum on smaller scales then the loop corrections caused by these quadratic 
terms become important. However, a consistent calculation of the halo correlation function or power spectrum at one 
loop level requires a model of the halo density field up to third order, including all the non-local terms [5j. This gives 
rise to several other terms that contribute to the power spectrum at the same order, each with a prefactor that may 
depend on the halo mass. A clean extraction of these non- linear bias coefficients from the power spectrum alone is 
very difficult and instead higher order correlations are needed to separate these terms. This paper is a first step in 
extracting the quadratic coupling terms from the bispectrum. The next step in this program is to extract the cubic 
order terms from the trispectrum: such an analysis will be presented elsewhere. For this reason a discussion of the 
implications on the two point function based on our results for the quadratic couplings cannot be complete. Still, 
there are some obvious implications of our results. One is that the non-local tidal tensor bias contributions must be 
included in the analysis of galaxy power spectrum and ignoring them may lead to incorrect conclusions. This effect 
is specially relevant for the broadband power, but will also have an effect on the position of the baryonic acoustic 
oscillations (BAO), just like the quadratic density term does if &2 ^ [35]. Our preliminary calculations suggest this 
effect is subdominant, with about a factor of 4 lower prefactor in front of b s 2 relative to the bi term, suggesting it will 
not strongly affect the measurement of the true BAO scale. 
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Appendix A: Matter density in LPT up to second order 

In the above derivation of the final halo density field arising from a local bias in Lagrangian space we used the 
Lagrangian displacement field and the second order matter density in SPT. In [37] it is shown that after expanding 
the exponential damping prefactor, the one loop power spectra in LPT and SPT are identical. This appendix shows 
that this equivalence is true also in terms of the fields if one expands the LPT expressions up to the desired order in 
the density field. Without these expansions, LPT was shown to contain a non-trivial resummation of SPT terms. 

For n-th order displacement field one has [37] 

(n) *M = -4 / §^f---J ^ {n) L{ Pl ,..., Pn )5(p 1 ^)...5{p n ,r 1 )(2K)HV)(k-Y J p), (Al) 
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with the kernels 



(2) 



L(Pl,P2) = 



3 Pl + P2 

7 |Pi +P2I 2 



Pl P2 

PlP2 



(A2) 



The density field in fc-space can be obtained upon Fourier transforming the Eulerian field in configuration space and 
using the Jacobian mapping [l + 15(03,77)] d 3 x — d 3 q 



5(k,rf) = J d 3 q exp [i fc • q] {cxp [i k ■ *(<?, 77)] - 1} 
We can now expand the exponential up to second order in the displacement field 
<5(fc,7?)« J d 3 g exp [ifc-q]ji 77)- i(fc • t;)) 2 J 

1 f d 3 fc' k-k'k- (fe-fe')(i) 



=ifc-( 1 )*(fc,77)+ife-( 2 )*(fc,?7) + i J 



(2tt) 3 (fc') 2 \k - k'\ 



S(k',r])^S(k-k',r]) 



-(i) 



<5(fc,77) 



(A3) 

(A4) 
(A5) 
(A6) 



Thus, expanding the exponential we see that LPT and SPT agree at second order in the density field. 
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